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ABSTRACT 



A data-driven method and system are provided for predict- 
ing limit cycle oscillations such as chatter to obtain a 
classifier signal which, in turn, may be utilized by a control 
method and system. The method and system utilize algo- 
rithms based on random field theory which are applied: to 
'Vibration" data to predict when chatter phenomenon may 
emerge during an ongoing machining process. Using real 
data from a turning operation on a lathe, the onset of chatter 
can be predicted so that sufficient time is available to deploy 
an automatic control strategy such as reduction or random 
variation of spindle speed to quench the "chatter". The 
method and system preferably use instantaneous scale of 
fluctuation. 
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METHOD AND SYSTEM FOR PREDICTING 
LIMIT CYCLE OSCILLATIONS AND 
CONTROL METHOD AND SYSTEM 
UTILIZING SAME 

TECHNICAL FIELD 

This invention relates to methods and systems for pre- 
dicting limit cycle oscillations and control methods and 
systems utilizing such predictions. 

BACKGROUND ART 

Limit cycle oscillations arise in many natural systems 
such as machine tool, brain electrical signals, return on 
stock, electrocardiogram, earthquake, vibrations of 
mechanical structures (such as bridges., airplane wings, 
automobile engines and suspensions among others) and 
communication networks, for example 

Chatter, a type of limit cycle oscillation, is a self-excited 
relative vibration between the workpiece and the cutting tool 
in common machining processes such as turning process on 
a lathe as described in Zhang, H., Ni. J. and Shi, H., 
"Machining Chatter Suppression by Means of Spindle Speed 
Variation", PROC. S. M. Wu Symposium, Vol. 1. pp. 
161-175. 1994. 

The development of charter causes a machining process to 
become unstable. This can result in poor surface finish, 
accuracy and reduced tool and other machine part life. 
Having to control machine tool chatter after it arises by 
traditional means such as by reducing spindle speed restricts 
the amount of material that can be removed in a turning 
process, for example. Therefore, development of chatter in 
a machining process reduces productivity and quality. 

The three major approaches to chatter prediction are 
described in Minis, I. E., Magrab, E. B. and Pandelidis, I. O., 
"Improved Methods for the Prediction of Chatter in Turning. 
Part 3: A Generalized Linear Theory". ASME Journal of 
Engineering for Industry, vol 112, pp: 2&-3S, 1990;TanseL 
L N. et al, Recognition of Chatter With Normal Networks," 
Int. J. Mack. Tools. Manufact., vol. 31. No. 4. pp. 539-552. 
1991; and the U.S. patent to Delio entitled "Method of 
Controlling Chatter in a Machine Tool," U.S. Pat No. 
5.170358. 

The presence and evolution of chatter can be monitored 
by measuring the vibration signals from the cutting tool 
using appropriately placed accelerometers. Such a signal is 
called the "chatter signal." The chatter signal can be ana- 
lyzed to discern the state of the machining process. A 
conceptual model of the chatter signal generation is shown 
in FTG. 1. As is clear, the chatter signal, X(t), will be a 
complicated random signal. 

The three approaches to manufacturing operations of 
relevance to chatter phenomenon can be categorized as (1) 
ftevent Approach (leads to conservative operational ratings 
and resultant low throughput); (2) Detect and Control 
Approach (leads to reworking or scrapping workpieces and 
resultant low quality and throughput); and (3) Predict and 
Prevent Approach (leads to optimal high operational ratings 
and resultant high quality and high throughput). 

SUMMARY OF THE INVENTION 

An object of the present invention is to provide a method 
and system for predicting limit cycle oscillations in a 
data-driven fashion. 

Another object of the present invention is to provide a 
method and system for predicting limit cycle oscillations in 
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a data-driven fashion to obtain a classifier signal which, in 
turn, is utilized in a control method and system of the present 
invention. 

Another object of the present invention is to provide a 

5 method and system for predicting limit cycle oscillations 
such as chatter wherein algorithms based on random field 
theory are applied to 'Vibration" data to predict when the 
chatter phenomenon may emerge during an ongoing 
machining operation. 

10 Still another object of the present invention is to provide 
a method and system for predicting limit cycle oscillations 
such as charter wherein algorithms based on random field 
theory are applied to 'Vibration" data to predict when the 
chatter phenomenon may emerge during an ongoing 

13 machining operation so far in advance that there is ample 
time to use such a prediction. Control strategies such as 
reduction or random variation of spindle speed to "quench" 
the chatter may then be employed in a control method and 
system of the present invention. 

20 Yet still another object of the present invention is to 
provide a method and system for predicting limit cycle 
oscillations in a machining operation to increase throughput 
and quality of machined workpieces. 

25 In carrying out the above objects and other objects of the 
present invention, a method is provided for predicting limit 
cycle oscillations. The method includes the steps of gener- 
ating a signal based on oscillatory behavior, statistically 
processing the signal to obtain a corresponding envelope 

30 signal and monitoring instantaneous changes in a function 
which is based on the envelope signal. The method also 
includes the step of generating a classifier signal when the 
instantaneous changes are less than a predetermined thresh- 
old value for a predetermined time period whereby the 

35 classifier signal is a prediction of the limit cycle oscillations. 
Still further in carrying out the above objects and other 
objects of the present invention, a control method is also 
provided utilizing all of the above-noted steps and further 
including the step of generating an output control signal 

40 based on the classifier signal 

In one embodiment, the function is an instantaneous scale 
of fluctuation which tracks changes in a scale of fluctuation 
over time. The instantaneous scale of fluctuation is a func- 
tion of an estimate of the energy distribution of the envelope 

45 signal and an inverse Fourier transform of the energy 
distribution of the envelope signal. 

In another embodiment the instantaneous scale of fluc- 
tuation is estimated also based on the energy distribution of 
the envelope signal and an estimate of the variance of the 

50 envelope signal. 

In yet still another embodiment, a simple approximation 
of the method of monitoring the instantaneous scale of 
fluctuation is provided by monitoring an instantaneous auto- 

55 correlation function or its approximate estimates and its first 
zero-crossing lag. 

Still further in carrying out the above objects and other 
objects of the present invention, systems are provided for 
carrying out the above-noted method steps. 

60 The ability to predict chatter allows for simultaneously 
maximizing production throughput surface finish of the 
workpiece and machine tool life in the industrial practice of 
machining process. The method and system of the present 
invention (also referred to as Limit Cycle Prediction Using 

65 Scale of Fluctuation (LP-SOF)) permits the prediction ; of 
chatter development over 300 milliseconds in advance, 
thereby providing an opportunity within the machining 
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process to avoid chatter and the concomitant reduction of 
productivity and quality. Therefore, the invention described 
herein should find practical application in any machining 
process where chatter phenomenon can compromise quality 
and productivity. s 

Chatter prediction using the method and system of die 
present invention relies on a purely data-driven approach 
and does not require the knowledge of various system model 
parameters a-priori. Newly-developed algorithms based on 
random field theory are applied to 'Vibration" data to predict 10 
when the chatter phenomenon may emerge during an 
on-going machining process. 

The method and system of the present invention makes 
the Predict & Prevent approach of FIG. 2 possible in the 
manufacture of high value-added components. FIG. 3 illus- 
trates the content of the module of FIG. 2. 

Using real data collected during machining operations, it 
has been shown that the onset of chatter can be predicted as 
much as 300 milliseconds in advance. This is ample time to 
deploy simple automatic control strategies such as reduction 
or random variation of the spindle speed to prevent the 
chatter (rather than "quench" chatter after it has fully 
developed, which can be much more difficult). The spindle 
speed control output signal, sc(n). of FIGS. 2 and 3 provide 
such control. 

Some of the commercial features of the method and 
system of the present invention are as follows: 
After some initial developmental work, the method and/or 
system of the present invention can be part of any new 
CNC machines or "retrofitted" on many existing 30 
machines. 

The potential for such a method and system in a high 
speed, high quality production system is obvious (ie.. 
increased throughput and quality). The financial impact 
of such an advance in the manufacture of high-value- 35 
added components such as pistons, computer disk drive 
parts, etc.. can be significant. 
The ability to predict chatter permits a certain extra 
"margin of safety" in the design and manufacturing 
stages (for example, the quality of raw material for 40 
manufacturing can be more viable) of machine tools. 
For the same level of design and manufacturing quality, 
the machine tool can be operated above its usual 
conservative rating. 
The method and system of the present invention provides 45 
machine tool chatter prediction in CNC machining systems 
utilizing newly developed methods in random field theory 
and Kalman-time frequency distributions. The machine tool 
chatter prediction is preferably accomplished by the specific 
method called LP-SOF. The ability to predict chatter allows so 
for simultaneously maximizing production throughput, sur- 
face finish of the workpiece and machine tool life in the 
industrial practice of machining process. 

The LP-SOF method uses the statistical signal processing 
methods of random field theory of local averages and 55 
Kalman time frequency distribution instead of the traditional 
methods of chatter prediction as described inTansel, L N. ct 
al, "Recognition of Chatter With Normal Networks," Int. J. 
Mack. Tools Manufact., Vol. 31. No. 4. pp. 539-552, 1991. 
to extract information from y(n), which is sampled version 60 
of X(t). The LP-SOF method can predict oncoming chatter 
over 300 milliseconds in advance. 

The above objects and other objects, features, and advan- 
tages of the present invention are readily apparent from the 
following detailed description of the best mode for carrying 65 
out the invention when taken in connection with the accom- 
panying drawings. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

FIG. 1 is a single block illustrating chatter signal genera- 
tion of the prior art; 

FIG. 2 is a single block illustrating output control signal 
generation of the method and system of the present inven- 
tion wherein Mode A is the chatter- free operation for a given 
feed rate and spindle speed and Mode B is for maximum 
spindle speed with chatter-free operation for a given feed 
rate; 

FIG. 3 is a schematic block diagram of the contents of 
Chatter P&P module of FIG. 2 wherein switch settings 
provide the Mode A or B operation; 

FIG. 4 is a schematic block diagram of the contents of the 
LP-SOF module or block of FIG. 3; 

FIG. 5 is a block diagram Sow chart illustrating a software 
implementation of the method and system of the present 
invention; 

FIG. 6 is a block diagram flow chart illustrating a hard- 
ware implementation of the method and system of the 
present invention; 

FIG. 7 are graphs of test data during normal operation (i:e. 
no chatter condition); 

FIG. 8 are graphs (similar to the graphs of FIG. 7 but with 
a different chatter signal amplitude scale) during various 
chatter conditions: (1) pre-chatter from 0 to 250 msec; (2) 
chatter onset from 250 to 450 msec; and (3) fully developed 
chatter from 450 msec to 1 sec; 

FIG. 9 are graphs which illustrate typical realizations of 
time series, AR A (2) with 6=17.9 and AR 2 (2) with 0=2; 

FIG. 10a is a graph with a 0=1 contour in a solid line and 
a circle, X 2 +y 2 in a dotted line; 

FIG. 10b is a graph with constant 6 contours; the two 
AR(2) random processes. A and B. are marked on the 6=3 
contour; 

FIG. 11 are graphs which illustrate typical realizations of 
time series, AR^(2) and AR^(2). both with 6=3; 

FIG. 12 is a graph with constant 6 contours; AR X and AR 3 
on 6=7 contour and AR 2 on 6=2 contour; 

FIGS. 13c, \3b and 13c are graphs which illustrate the 
results of simulation studies involving 3 AR(2) process, 
respectively. AR A with 9 A =7, AR 2 with 63=2, and AR 3 with 
63=7 ; typical realizations are shown in the top row and 
time-varying scale of fluctuation estimates in the bottom 
row; mean values of T(o). shown as dotted horizontal lines, 
are 7.2. 2.2 and 6.6 for AR^ AR 2 and AR 3 . respectively; and 

FIGS. 14a, 14b and 14c are graphs which illustrate 
eigenvalues of the averaged normalized autocovariance 
matrices of AR X . AR 2 and AR 3 . respectively. 

BEST MODE FOR CARRYING OUT THE 
INVENTION 

The method and system of the present invention can be 
best explained by defining precisely the functions of each of 
the blocks in FIG. 4 and specifying the proper choice of 
parameters therein. FIG. 4 is a detailed block diagram of die 
components of the LP-SOF module of FIG. 3. 
Signal Conversion 

A sampled chatter signal. y(n). is converted to its "enve- 
lope signal", e(n) at block 10. The y(n) signal's analytic 
signal counterpart is obtained by Gabor's construction and 
the magnitude of the analytic signal is die envelope. e(n). 

Analytic signal. a(n)=y(n>fj HT[y(n)), where M j" is the 
V(-l) and HT(] is the Hilbert transform operation. Then, the 
envelope. e(n)=Ja(n)l where II is the magnitude of the 
complex number, a(n). The envelope. e(n) is real and low 
pass. 
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Kalman-TFD around 7. Thus T(n) is relatively insensitive to operat- 

As described in Madhavan. P. G, and Williams, W. J.. "A ing conditions during normal operation and fully devel- 

TLme Varying Model for the Estimation of Time-Frequency oped chatter. 

Distribution of Signals in the Presence of Noise". fEEE c) In FIG. 8, the value of T(n) is significantly and 

Signal Processing Letters, 1996 (submitted), a time- varying. 5 consistently below 5 during pre -chatter condition. It 

stochastic, state space model of e(n) is used to estimate die can be concluded that the value of T(n) staying below 

energy distribution of e(n) in the time-frequency plane at 5 is a predictor of chatter mat develops fully after 

block 12: approximately 400 msec 

Therefore, this working model demonstrates that if the 

X(n + i)= io chatter signal is processed as defined in FIG. 4, the value of 

r T(n) can be monitored to observe if it falls below 5. If T(n) 

, \ t j \ a. i T Mt «m-i> 1 stays below 5 for over 10 milliseconds, predict fully devel- 

M<«) fc- m Li...^ " J oped chatterto emerge in over 100 milliseconds. 

FIG. 5 is a block diagram flow chart of the steps pa- 
It has been shown that X^(n) can be estimated in the 15 formed by a properly programmed microcomputer or con- 
presence of noise using the Kalman filter algorithm and that troller to implement the method and system of the present 
SQui) is a good estimate of the instantaneous TFD. In the invention. The steps correspond to the blocks of FIG. 4 and 
present example. M=10. the corresponding description thereof as noted above. 
Scale of Fluctuation (i.e.. block 14) FIG. 6 is a schematic block diagram of a hardware 
As described in Vanmarcke. E., "Random Fields: Analysis 20 implementation of the method and system of the present 
and Synmesis".A//TPwjj; Cambridge. Mass. 1988. scale of invention. The blocks of FIG. 6 also correspond to the 
fluctuation, 0. captures some important statistical properties blocks of FIG. 4 and the corresponding description thereof 
of time series such as "characteristic or correlation length." as noted above. 
The scale of fluctuation, 9. is defined mathematically as Alternative Embodiments 

follows: 25 An alternate method of estimating the instantaneous scale 



of fluctuation. T(n). is as follows: 



I S(k.n) 

where G(0) is the zero- frequency value of the one-sided 30 fel 
power spectral density of e(n) and B(0) is the zero-lag value 

of the autocorrelation function of e(n). where S(0ji) is the zero-frequency slice of the Kalman-TFD 

Anew measure of the present invention called the "instan- an d the denominator is an estimate of the variance of die 

taneous scale of fluctuation" "tracks" the changes of 6 over signal. Even though T^n) does not have the usual meaning 

time. Based on the Kalman-TFD. a new estimate of the 35 0 f scale of fluctuation. T^fn) will be normalized such that it 

instantaneous scale of fluctuation is defined as follows: always lies between 0 and 1. This may be advantageous in 

certain applications where threshold monitoring is utilized. 

7fri)= S(0,n) As a crude and simple approximation of the LP-SOF 

method of monitoring the instantaneous scale of fluctuation. 

40 the instantaneous autocorrelation function. r(0.n). (or its 

where S(0.n) is the zero-frequency slice of the Kalman-TFD approximate estimates such as short-time windowed auto- 

and r(Oa) is the zero-lag slice of instantaneous autocorre- correlation function) can be utilized. The first zero-crossing 

lation function, r{m,n)=IFT|S(kji)] where IFT[] is the ^ T<> 0 f the instantaneous autocorrelation function. r(0,n), 

inverse Fourier transform (Le.. block 16). can be monitored and significant reduction of its value can 

Threshold Detection 45 potentially be used as a means of predicting the onset of 

Referring again to FIG. 3, the estimated instantaneous limit cycle oscillations, 

scale of fluctuation, T(n), is continuously monitored. If T(n) Ramifications 

is less than 5 far more than 10 milliseconds, fully developed The LP-SOF method can be used to predict the appear- 

chatter is predicted in 100 milliseconds or later. This means a nce or disappearance of limit cycle oscillations in many 

that the CNC machine control system will have 100 milli- 50 natural systems as noted in the Background Art of the 

seconds or more to institute chatter control measures. present application. 

Working Models Scale of Fluctuation — An Invariant Property of a Class: of 

The LP-SOF method was tested with real data from Random Processes 
turning process on a lathe during normal operation and Vanmarcke developed a new theory of random fields 
chatter development Results are shown in FIGS. 7 and 8. 55 b asc( j on foe concept of local averaging and the character- 
Both Figures show the chatter signal, y(n) on the top panel ization of the second-order properties of the random fields 
and the smoothed estimate of instantaneous scale of using the variance function, Certain asymptotic properties of 
fluctuation. T(n). in the bottom panel. the variance function lead to the definition of a scalar called 

The following features are important to note in FIGS. 7 the "scale of fluctuation. 9". which has many interesting 

and 8: 60 properties. A newly-discovered property is the invariance of 

a) The amplitude scales for chatter signal are different; the e for certain class of random processes. This feature is 
amplitude of x(n) in FIG. 7 during normal operation is investigated for second-order auto-regressive random pro- 
of the same order as the amplitude of prc-chatter signal cesses and it is shown that constant -8 contours for AR(2) 
in the interval from 0 to 250 msec in FIG. 8. processes are approximately circles with centers on die 

b) In FIG. 7, the value of T(n) stays relatively stable 65 non-negative real axis with radius less than or equal to half, 
around 7. In FIG. 8. when chatter is fully developed A non-parametric method of estimating time-varying 8 is 
after 450 msec. T(n) value is again approximately developed using the time-varying model-based time- 
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frequency distribution. Simulated machine tool chatter data 
show that 9 captures a very important invariant property 
related to the unit-standard deviation contour volume of the 
joint probability density function of the chatter process. 
Tests with real vibration data from machine tools before and 
during chatter show that estimated 6 may permit on-line 
prediction of chatter development many hundreds of milli- 
seconds in advance. 
I Introduction 

Vanmarcke [11] introduced a comprehensive theory of 
random fields which extends very elegantly to multi- 
dimensional cases. His 'random field theory of local aver- 
ages** captures the effect that local averaging has on a 
homogeneous random field. The quantification of the effects 
of local averaging leads to a function which characterizes 
the second-order properties of the random field called die 
'Variance function." A scalar called the "scale of fluctuation** 
derived from the variance function has many nice 
interpretations, some historic and some new. The scale of 
fluctuation can be considered to be similar to other scalars 
derived from multidimensional probability density functions 
such as correlation or Shannon entropy. The information that 
the scale of fluctuation provides is different from correlation 
or entropy; in the case of a time series, the scale of 
fluctuation is a measure of the 'time scale** of a random 25 
process over which the correlation structure of that random 
process is characterized [9. 11]. 

After reviewing the basic concepts of random field theory 
and scale of fluctuation, this paper introduces a new invari- 
ance property of scale of fluctuation. It is shown that in the 30 
case of Gaussian random process, this invariance property 
captures the 'Volume" of the joint probability density func- 
tion which may have significant practical applications. A 
simulated data example is discussed later to develop the 
ideas behind the invariance property. An important notion of 35 
time-varying scale of fluctuation is introduced and a new 
estimation method using Kalman filtering is developed. The 
paper concludes with preliminary results from a practical 
application of time-varying scale of fluctuation for the 
prediction of machine-tool chatter. Tests with simulated 40 
machine-tool chatter data provide the possible rationale for 
the use of time-varying scale of fluctuation for the prediction 
of machine-tool chatter. 
H. Theory Review 

Consider a zero-mean, stationary, real random process, 
x(t) with a variance of a 2 and the following definitions of its 
second-order properties: 

1) Correlation Function: 



tf(T)=£lWm)J 



2) Normalized Correlation Function: 



3) Spectral Density Function (s.dX): 

5(<D)= ~5T | 

4) One-sided s.d.f.; 

G(cd)=25(cd); u>20 



50 



55 



60 



8 



5) Normalized s.dJ. (unit area): 



The concept of variance function is developed by con- 
sidering the following moving average process. Xj(t): 



For larger averaging window, T, the variance, a x 2 , of the 
corresponding moving average process will be smaller than 
O" 2 . The variance function is defined as follows: 



Vanmarcke showed that as T increases, TfT). has a charac- 
teristic shape and that the limiting value of Ty(T) is very 
meaningful. In fact, the limit is called the "scale of 
fluctuation. 6** and the limit exists [11] if the slope of the 
normalized spectral density function at £0=0 is zero. 



■j: 



p<T>fr = n*(0) 



From equation (1), 6 is die "area** under the normalized 
correlation function or proportional to the zero-frequency 
ordinate of the normalized spectral density function. Many 
interpretations for G exists and is variously known as ''char- 
acteristic length**, "correlation length** or "duration of per- 
sistence of trends.** For a real-valued autoregressive random 
process. X(t), or order P, using results developed by Van- 
marcke [11], the scale of fluctuation, 6^ can be obtained as 
follows: 



9, = 6* 



45 



a? 



<2) 



where !//(0)P = 



[ 4« 1 



Here, H(co) is the transfer function of the linear-time invari- 
ant system relating white noise process to the AR(P) process, 
6* is the scale of fluctuation of the white noise process 
(0^=1), o w 2 is the variance of the white noise process 
(o„ 2 =l), a, 2 is the variance of the AR(P) process and {a,}. 
i=l to P are the AR coefficients. 

To develop a better understanding of 9 consider two 
discrete-time AR(2) processes with model equations as 
follows-<l) AR A (2): x^l^n-lH^ x^n-^Hw, 
(n) and (2) AR 2 (2): x 2 (n)==0.6x 2 (n-lH).2x 2 (n-2>+w 2 (n). 
Topical 100-point data sequences are shown in FIG. 9. Using 
equation (2). the scale of fluctuation can be calculated 
knowing that [4], 



65 



1 +aj 
1 



1 



P) 



where a t . aj are the AR coefficients {-1.8. 0.82} and {-0.6. 
0.2}. For AR J (2). 6^17.9 and for AR 2 (2), 82=2. The intcr- 
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pretations of 9 as "correlation length" or "duration of 

persistence of trends" are clear by considering the plots in x^i^n^n); WPUM . - - WO . . . *Mf M 

FIG. 9. The ARi(2) data is seen to be a narrowband signal 

with correlation function that will persist for a long time X«) = fc£*(n) + «(«); ^ 

(9,= 17.9) whereas the AR 2 (2) data has a correlation tunc- 5 

tion that dies out rapidly (6 2 =2). r 
m Invariance Property & = I? I i ...V"^" *...e~^~ ,,( *~ ,, ] 

For a real-valued AR(2) process, consider the expression 
for 6. If the pole positions corresponding to real coefficient 10 For arbitrary choice of the matrix. A„, and the vector <fr„. 
{a 4 , aj are (x±jy). the expression for 6 can be written as cquati ons (5a), (5b) is sometimes referred to as a time- 
follows: variable parameter (TVP) model of a stochastic time series 



(l+j«» + y»)[(l-Z»+jri + ^] 15 



[13]. In the noise-free observation case. y(n) can be written 

as: 

(6) 



The "9-1" constant-0 contour is plotted in FIG. 10a. It = $&(«)= -j^ ^ c " X M 

shows the plot in the first quadrant of the z-plane with a 

serai-circle for comparison. By setting 9=1 in equation (4), Equation (6) can be seen to be the Goertzel algorithm [8] 
we obtain the quartic equation. (x J +y 2 ) 2 +x 2 V-2x=0. M where XJWl) is the kth coefficient of ite M-^l DIT die 

shown in FIG. 10a. For comparison, the circle shown in FIG. fijute-duration sequence {y(n) ; 0£nS(M-l)} and X(n) is 

10a has the quadratic equation, <x-fc) 2 +yMV*) 2 . If we the vector, elements of which are the M-point DFT 

consider the 8=1 contour approximately circular, in FIG. coefficients. **<M). K-i to M. 

;T : " y y . A ; - ' An alternate way to consider equation (6) is by compan- 
ion we can see that constant-e contours for 0>1 are son ^ ^ ^ presentation when y(n) is 
approximate circles of smaller diameter and with centers * restricted l0 M widc _ SC osc stationary [10]. The Cramer 
along the real axis closer to the unit circle. Any pole pairs represe ntation of y(o) can be written as: 
that lie on the same constant-6 contour will have different 
{a x . a 2 } coefficients but the same 9. For example, consider 

the two AR(2) processes marked on the 9=3 contour in FIG. ^ _ 1 f n 

10£> (only the pole in the first quadrant is marked). The first 211 J -n 

process, AR^(2), has poles at (0.9±j0.26) and the second 

process, AR^). has poles at (0.35±j0.17). The correspond- nc increment process {dz(co)} has the properties that its 

ing AR coefficients are {-1.8. 0.8776} and {-0.7. 0.154}. energy at different frequencies is uncorrelated and that the 
respectively. Typical 100-point data sequences are shown in 35 expected value of ld2(co)l 2 defines the spectrum S y (co)dca 

FIG. 11. As is clear, AR A (2) and AR^2) are very different Analogous to ldZ(co)l 2 , for the nonstationary case in equation 

processes with different correlation functions (or spectral (2), we can define a time-varying spectrum based on 

densities), yet the same 9. The significance of this invariance X(n). We proceed as follows. 

property will become clearer when we consider its practical In equation (5), assume that e(n) and w.(n) are zero-mean 
application in the case of machine-tool chatter signal analy- 40 Gaussian white noise sequences such that E|e(k)w_(l)]=Q for 

sis and related simulation studies. all k and 1 and the initial "state" X(0) is independent of e(k) 

At this time, it is not clear what invariant physical £^ 
property of Jhe random process is being captured by the J observation g y( ^ me 

scale of fluctuation. However, it is interesting to consider ^^mp** estimates X(n) of the DFT coefficients at 
further the linear time invariant system in the modeling of variations of the DFT coefficients over time can 

AR processes. In this simple, 2-order linear system case, we be used to estimatc me time-frequency distribution [ 1. 121. 

can compare the constant-9 contour to other traditional We Mnt me time-varying model based time-frequency 

contours such as constant-^ (damping ratio) and constant-o) n distribute,, (TVM-TFD) as : 
(natural frequency) contours [3]. The constant-*; contours 
are logarithmic spirals and the constanMu„ contours are at 50 

right angles to the logarithmic spirals, both distinctly dif- S(t/i)=Uf»P for/i,t=U, . . . N (7) 

ferent from the constant-9 contour Our conjecture is that all ^ t^.^ estimatcs ^ ^ spec tral 

Gaussian random processes that be on the constant-9 con- dcn$ . ^ ^ u k ^ ^ m By 

tour have similar unit-standard deviation contour volume txicn 6isig equation (1) by the use of time-varying norraal- 

This possibility will be explored further for the case of . m& density ((Dt) a time.varying 9 can be 

simulated machine-tool chatter signal analysis. dcfincd M ^ WQ?olXioti ^ to me ZC ro-frequency ordinate 

IV. Time- Varying "6" of Discrete-Time Signals 0 f the normalized time-varying spectral density function. 

In recent years, significant advances have occurred in the i.e.. 9(t>=7Cg(0.t). For the discrete-tie case, based on the 
field time-varying spectra or 'time-frequency distributions" ^ TVM-TFD. we define an estimate of the time-varying scale 

fl). The time-frequency distributions (TFDs) allow you to of fluctuation. T(n), as follows: 
estimate the power spectral density as it varies over time, 
primarily for deterministic signals. These concepts have S(Q/i) 
been extended to random signals by the de velop ment of 
time-varying model based TFD called "TVM-TFD" (6|. 65 Here, r(ra.n) is die instantaneous autocorrelation sequence 

Consider a general state space model of a time-varying and r(0ai) is the instantaneous mean-squared value of the 

stochastic discrete-time signal. y(n), as given below: discrete-time random process, y(n). 



(8) 
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V. Time- Varying <4 6" — Simulation Studies 
The properties of the estimate of the time-varying scale of 
fluctuation. T(n). will closely resemble the properties of the 
TVM-TFD [61. The use of Kalman filtering for TVM-TFD 
assures optimality of estimates in the mean-square sense 
which translates into good noise performance. The good 
localization and cross-term properties of the TVM-TFD 
results in S(O.n) being a good time-varying estimate of the 
zero-frequency ordinate of the time-varying spectral density 
function. In this section, we demonstrate some of these 
properties through simulation studies. H\c simulation data 
are chosen with the additional aim of providing a basis for 
die explanation of real data tests to be undertaken in Section 
VL 

The pole positions of the three AR(2) processes are shown 
in FIG. 12. The processes are (1) ARj, with poles at 
(0.6±j0.095) and corresponding AR coefficients {-1.2, 
0.369}. (2) AR 2 . with poles at (0.6±j0.46) and corresponding 
AR coefficients {-1.2, 0.5716} and (3) AR 3 , with poles at 
(0.95±j0.15) and corresponding AR coefficients {-1.9, 
0.9221}. The G of the three processes are calculated by 
equation (2) and marked on the constant -6 contours in FIG. 
12. The 6 for AR,. and AR 3 . G 3 are equal (G,=0 3 =7) 
whereas for AR 2 , 02=2. The FIGS. 13a, 136 and 13c show 
a realization of 1000 points each of the three time series as 
well as plots of the estimates of the time-varying scale , of 
fluctuation* T(n), for each time series. By visual inspection 
of the time series plots, it can be seen that the time-domain 
(amplitudes) and frequency-domain (broad-band nature) 
features of A, and AR 2 are very similar and AR 3 is most 
dissimilar whereas T(n) of AR X and AR 3 are similar 
(estimate close to the theoretical value of 7) and AR 2 is most 
dissimilar (estimate close to the theoretical value of 2). It can 
be observed that T(n) "tracks" the theoretical 6 value for the 
three processes quite closely. The mean values of T(n) 
(shown by dotted horizontal line) in FIGS. 13a, 13b and 13c 
are 7.2. 2.2 and 6.6, showing good agreement with the 
theoretical 9 values. 

Comparing FIGS. 13a and 13c, it is striking to note that 
whereas their G values are equal (=7), the appearance of the 
time series is quite different The AR 3 series is over 5 times 
in amplitude compared to AR r The AR 3 series is also quite 
band-limited compared to the broad-band nature of AR A . 
Note that these features are entirely consistent with the pole 
positions shown in FIG. 12 for these two processes. 
Conversely, from FIGS. 13a and 13b, it is striking to note 
that whereas their G values are different (7 and 2, 
respectively), the appearance of the time series is quite 
similar. This leads to the conclusion that 8 captures some 
property of the random process different from its first or 
second order statistical properties. In other words, even 
though the definition of 6 is based on the autocorrelation or 
spectral density function* G may capture one (or some small 
set of) specific feature of the joint probability density 
function. The intuition related to the real-data application to 
be discussed in Section VI leads to the possibility that the 
specific feature 8 captures may be related to the volume of 
the unit standard deviation contour of the joint probability 
density function. We explore this possibility as follows. 

Consider an n-dimensional random vector, y, representing 
a Gaussian discrete-time random process. The joint prob- 
ability density function of y is given by the following 
equation: 
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to) = ^ «p [ - t O " " ] 

5 Here, u is the mean vector and C is the autocovariance 
matrix. For the zero-mean case, the unit standard deviation 
(u.s.d.) contour is given by the equation. y r C-ly+l [see 2. 
for example]. In the 2-dimensional case, the equation is that 
of an ellipse with the semi major and minor axes equal to! ? 

10 times the product of the square roots of eigenvalues of C the 
autocovariance matrix [2]. Generalizing the n-dimensional 
case, the volume. V, of the unit standard deviation contour 
(a hyper-ellipse) is given by the following expression, where 
Xf is an eigenvalue of C. 

15 



In our simulation studies, the u.s.d. contour volumes were 
estimated for the three random processes from each of their 
(100x100) averaged, normalized autocovariance matrices. 
In FIGS. 14(a). 14(b) and 14(c). the largest 50 eigenvalues 
in descending order are plotted for AR t . AR 2 and AR 3 . It can 
be seen that the number of significant eigenvalues are 
different for each process. The selection of the number of 
significant eigenvalues is generally problematic. To aid in 
the selection process, we proceeded as follows. For AR i and 
AR 3 in FIGS. 14a and 14c, we calculated the "change of 
slope" of the ordered eigenvalues to aid in the detection of 
the "knee" in the eigenvalue profiles. The largest index for 
which the change of slope was grater than 1 was judged to 
be the location of the knee. The number of significant 
eigenvalues (index number. K) was chosen as one index 
number less than the knee location. Such a procedure 
yielded values of 7 and 6 for AR, and AR 3 . For AR 2 in FIG. 
14b, a "knee" is not detectable; K was chosen as 25 at which 
a sudden change of the magnitude of the eigenvalue can be 
detected. The u.s.d. contour volume of each joint probability 
density function is calculated by equation (10) where n is 
replaced by K for that random process. 



TABLE 1 



45 



Comparison of Unit Standard Deviation Nfolume 


AR Process 


0 


Number of Significant 
Eigenvalues, K 


Unit Standard Deviation 
Contour Vohnne 


AR, 


7 


7 


2.0 x 10 3 


AR, 


2 


25 


2.03 x 10* 


AR, 


7 


6 


6.4 x 10 3 



50 



The comparison in Table 1 shows that ARj and AR 3 have 
similar 0, "degrees of freedom" and volumes whereas AR 2 
has a volume 3 orders of magnitude larger. This observation 
justifies the conjecture mat 0 may capture a measure related 

55 to the volume of the unit standard deviation contour. The fact 
that small 0 implies large volume in which the random 
process exists suggest an intuitively appealing explanation 
of the results of the real-data application to be discussed in 
the next section. 

60 VT. An Application Of Time-Varying "6 M : Prediction Of 
Machine Tool Chatter 

Machine tool chatter is a self-excited relative vibration 
between the workpiece and the cutting tool in common 
machining processes such as turning process on a lathe [7]. 

65 The presence and evolution of chatter can be monitored by 
measuring the vibration signals from the cutting tool using 
appropriately placed accelerometers. Such a signal is called 
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the "chatter signal." The chatter signal can be analyzed to 
discern the stale of the machining process and this informa- 
tion can be used to predict the possible future development 
of chatter. 

In FIGS. 7 and 8. the chatter signal and the corresponding 
estimates of the time-varying scale of fluctuation, T(n), are 
shown. In FIG. 7, the condition of the machine is "normal 
operation" It can be noticed that the time series and T(n) are 
comparable to AR, in the simulation studies in Section V 
(compare FIGS. 7 and 13a). In FIG. 8. the first 300 milli- 
seconds or so is the **pre-chatter" period and the rest is the 
•fully -developed chatter" period. The "p^e-chatte^' , period is 
comparable to FIG. 13fe and 'fully-developed chatter." to 
FIG. 13c. 

The following features are important to note in FIGS. 7 
and 8. The amplitude scales for chatter signal are different; 
the amplitude in FIG. 7 during normal operation is of the 
same order as the amplitude of pre-chatter signal in the 
interval from 0 to 250 msec in FIG. 8. In FIG. 7. the value 
of T(n) stays relatively stable around 7. In FIG. 8, when 
chatter is fully developed after 450 msec, T(n) value is again 
approximately around 7. Thus T(n) is relatively insensitive* 
to operating conditions during normal operation and fully 
developed chatter. In FIG. 7, the value of T(n) is signifi- 
cantly and consistently below 5 during pre-chatter condition. 
It can be concluded that the value of T(n) staying below 5 
is a predictor of chatter that develops fully after approxi- 
mately 400 msec. 

The reduction of the value of T(n) during pre-chatter 
period as a predictor of future chatter development can be 
explained in the light of observations made in Section V. The 
random process during pre-chatter condition behaves similar 
to AR 2 in Section V, which as shown in Table 1, goes 
through an increase in "degrees-of-frcedom" or its u.s.<L 
contour volume. This expansion of the machine system's 
phase-space makes it vulnerable to be trapped into certain 
fundamental modes of the system which will result in energy 
build up in a narrow frequency range and resultant nonlinear 
limit cycle type of oscillations characteristic of chatter. 
VTL Conclusions 

The invariance property of scale fluctuation. 6, for certain 
class of random processes was introduced For second-order 
auto-regressive random processes, it was shown that 
con stan t-6 contours are approximately circles with centers 
on the non-negative real axis with radius less than or equal 
to half. It should be noted that the invariance property of 6 
discussed in this article is distinct from a quantity Van- 
marcke has defined as "invariant-6" of a random process 
which is the product of the variance and the 0 of the random 
process 111). It can also be seen that in the second order 
system case, the contours of constant "invariant-6" is quite 
different from the approximately circular constant 0 con- 
tours discussed in this article. 

A Kalman filter based non -parametric method of estimat- 
ing time-varying 6 was developed using the time-varying 
model-based time-frequency distribution. Simulated data 
showed that 9 captures a very important invariant property 
related to the unit-standard deviation contour volume of the 
joint probability density function of the Gaussian random 
process. Tests with real vibration data from machine tools 
before and during chatter show that estimated 9 may permit 
on-line prediction of chatter development many hundreds of 
milliseconds in advance. 
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While the best mode for carrying out the invention has 
been described in detail, those familiar with the art to which 
this invention relates will recognize various alternative 
designs and embodiments for practicing the invention: as 
defined by the following claims. 
What is claimed is: 

1. A method for predicting limit cycle oscillations in a 
natural system, the method comprising the steps of: 

measuring oscillatory behavior in the natural system to 
obtain a first signal based on the oscillatory behavior; 
statistically processing the first signal to obtain a corre- 
sponding envelope signal; 
monitoring instantaneous changes in a function which is 

based on the envelope signal; and 
generating a classifier signal when the instantaneous 
changes are less than a predetermined threshold value 
for a predetermined time period whereby the classifier 
signal is a prediction of the limit cycle oscillations; 

2. The method as claimed in claim 1 wherein the limit 
60 cycle oscillations are undesirable vibrations. 

3. The method of claim 2 wherein the undesirable vibra- 
tions are machine tool chatter. 

.4. The method of claim 1 wherein the step of monitoring 
includes the step of estimating an instantaneous scale of 
65 fluctuation of the envelope signal and wherein the instanta- 
neous changes are based . on the instantaneous time- 
frequency distribution. 
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5. The method of claim 4 wherein the instantaneous 
time- frequency distribution is a Kalman time-frequency 
distribution. 

6. The method of claim 1 wherein the first signal is a 
vibration signal. 

7. The method of claim 1 wherein the first signal is a 
sampled vibration signal. 

8. The method of claim 1 wherein the function is an 
instantaneous scale of fluctuation. 

9. The method of claim 1 wherein the function is an 
autocorrelation function 

10. A method for controlling limit cycle oscillations of an 
object the method comprising the steps of: 

generating an electrical signal based on oscillatory behav- 
ior of the object; 

statistically processing the electrical signal to obtain a 
corresponding envelope signal; 

monitoring instantaneous changes in a function which is 
based on the envelope signal; 

generating a classifier signal when the instantaneous 
changes are less than a predetermined threshold value 
for a predetermined time period whereby the classifier 
signal is a prediction of the limit cycle oscillations; and 

controlling the limit cycle oscillations of the object based 
on the classifier signal. 

11. The method as claimed in claim 10 wherein the step 
of controlling includes the step of generating an output 
control signal adapted to control a machine tool 

12. The method of claim 11 wherein the output control 
signal is a speed control signal. 

13. The method of claim 12 wherein the object is a 
machine tool having a spindle and wherein the speed control 
signal is a spindle speed control signal. 

14. A system for predicting limit cycle oscillations, the 
system comprising: 

a detector for generating a first signal based on oscillatory 
behavior, 

a signal processor for statistically processing the signal to 
obtain a corresponding envelope signal; 

a monitor for monitoring instantaneous changes in a 
function which is based on the envelope signal; and 

a signal generator for generating a classifier signal when 
the instantaneous changes are less than a predetermined 
threshold value for a predetermined time period 
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whereby the classifier signal is a prediction of the limit 
cycle oscillations. 
15. The system of claim 14 wherein the limit cycle 
oscillations are undesirable vibrations. 
5 16. The system of claim 15 wherein the undesirable 
vibrations are machine tool chatter. 

17. The system of claim 14 wherein the monitor includes 
an estimator for estimating instantaneous scale of fluctuation 
of the envelope signal and wherein the instantaneous 
changes are based on the instantaneous time-frequency 

1 distribution. 

18. The system of claim 17 wherein the instantaneous 
time-frequency distribution is a Kalman time-frequency 
distribution. 

19. The system of claim 14 wherein the first signal is a 
15 vibration signal. 

20. The system of claim 14 wherein the first signal is. a 
sampled vibration signal. 

21. The system of claim 14 wherein the function is an 
instantaneous scale of fluctuation. 

20 22 The system as claimed in claim 14 wherein the 
function is an autocorrelation function. 

23. A system for controlling limit cycle oscillations of an 
object, the system comprising: 

a detector for generating an electrical signal based on 
23 oscillatory behavior of the object; 

a signal processor for statistically processing the electrical 

signal to obtain a corresponding envelope signal; 
a monitor for monitoring instantaneous changes in a 
30 function which is based on the envelope signal; 

a first signal generator means for generating a classifier 
signal when the instantaneous changes are less than a 
predetermined threshold value for a predetermined 
time period whereby the classifier signal is a prediction 
35 of the limit cycle oscillations; and 

a second signal generator for generating an output control 
signal based on the classifier signal. 

24. The system as claimed in claim 23 wherein the output 
control signal is adapted to control a machine tool. 

40 25. The system as claimed in claim 24 wherein the output 
control signal is a speed control signal. 

26. The system as claimed in claim 25 wherein the object 
is a machine tool having a spindle and wherein the speed 
control signal is a spindle speed control signal. 
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